function ds = xeoms(t, s, mu, cexh)
if t < 1
    alpha = atan2(s(4), s(6));
else
    alpha = atan2(-s(4), -s(6));
end
beta = 0;

ds = zeros(7, 1);
ds(1) = s(4);
ds(2) = s(5)/s(1);
ds(3) = s(6)/(s(1)*sin(s(2)));

ds(4) = s(5)^2/s(1) + s(6)^2/s(1) - mu/s(1)^2 + s(7)*sin(alpha)*cos(beta);
ds(5) = -s(4)*s(5)/s(1) + s(6)^2/(s(1)*tan(s(2))) - s(7)*sin(beta);
ds(6) = -s(4)*s(6)/s(1) - s(5)*s(6)/(s(1)*tan(s(2))) + s(7)*cos(alpha)*cos(beta);

ds(7) = s(7)^2/cexh;
% ds(7) = 0;
end